***Code used to create the table and figures presented in the manuscript:
***Incentivizing COVID-19 vaccination among racial/ethnic minority adults in the United States: $209 per dose could convince the hesitant 
***Authors: Kevin Chen, Marta Wilson-Barthes, Jeffrey E. Harris, Omar Galárraga
***Latest analysis date: 12/06/2022; 01/14/2022

clear
capture log close
set more off

use "P:\Omar Projects\Incentives to Vaccinate Against COVID-19\Final Datasets - receibved 21Sep2021\I-VAC Incentives to Vaccinate Against Covid-19_Clean.dta"
log using "P:\Omar Projects\Incentives to Vaccinate Against COVID-19\Health Economics Review\HER - Resubmission due 4Nov2022\HER - R&&#2\IVAC Survey WTA_Ordered Logit Analysis.log", replace


***Create Table 1. Descriptive Statistics
		*use table1 command
		table1_mc, by(Vaccine) vars(UserLanguage cat \ Age cat \ Gender cat \ QUOTA cat \ AMT cat \ NEWAMT cat) total(before)

			*recode OtherPeople so that 3 = highest confidence (yes)
			gen Others2 = 3 if OtherPeople==1       //yes WTA amount
			replace Others2 = 2 if OtherPeople ==3  //unsure if WTA amount
			replace Others2 = 1 if OtherPeople ==2  //no not WTA amount

			*recode You so that 3 = highest confidence (yes)
			gen You2 = 3 if You==1       //yes
			replace You2 = 2 if You ==3  //unsure
			replace You2 = 1 if You ==2  //no

			*create 3 level variable for gender
			gen sex = . 
			replace sex = 1 if Gender == 2 |  Gender == 5 //male
			replace sex = 2 if Gender ==7 |  Gender == 6 //female
			replace sex = 3 if sex == . //all other gender categories 
			tab Gender sex

			*create numeric langauge var
			gen language = 1 if UserLanguage=="EN"
			replace language = 0 if UserLanguage=="ES"

		**table1 
		table1_mc, by(Vaccine) vars(UserLanguage cat \ Age cat \ sex cat \ QUOTA cat \ BLACKAFRICANAMERICAN cat \ HISPANIC cat \ ChooseNo cat) total(before)

		
***look at distributon of the random amount of the initial incentive
	**Randomly generated number (1 thru 10), determines compensation amounts presented
	tab AMT		//367 individuals divided into 10 groups, amounts shown for the question gauging OTHER PEOPLES' WTA
	tab NEWAMT 	//104 individuals divided into 10 groups, amounts shown for the question gauging RESPONDENT'S WTA
	tab Vaccine //263 already vaccinated 
	
***Figure 2A OtherPeople Undifferentiated by Vax - NO PERCENTAGES [KC's code]
gen amt_cat_2 = recode(AMT, 10, 20, 30, 40, 50)

label define cat_bins 10 "0-10" 20 "11-20" 30 "21-30" 40 "31-40" 50 "41-50"
label values amt_cat_2 cat_bins

set scheme s1mono
graph bar (count) [fw=amt_cat_2], over(OtherPeople) over(amt_cat_2) percent ///
stack asyvars bar (3, bcolor(gs14)) ///
bar (2, bcolor(gs10)) bar (1, bcolor(gs6)) blabel(bar, position(inside) format(%9.1f) color(black)) ///
yscale(r(-5 100)) ytitle("Proportion", size(medium)) ylabel(, angle(h)) ///
b1title("Incentive Amounts in USD") ///
legend(title("Response to Willingness-to-Accept", size(small))) ///
name(bar_other, replace)
	
***Figure 2B You Undifferentiated by Vax - NO PERCENTAGES [KC's code]
set scheme s1mono
graph bar (count) [fw=amt_cat_2], over(You2) over(amt_cat_2) percent ///
stack asyvars bar (3, bcolor(gs14)) ///
bar (2, bcolor(gs10)) bar (1, bcolor(gs6)) blabel(bar, position(inside) format(%9.1f) color(black)) ///
yscale(r(-5 100)) ytitle("Proportion", size(medium)) ylabel(, angle(h)) ///
b1title("Incentive Amounts in USD") ///
legend(title("Response to Willingness-to-Accept", size(small))) ///
name(bar_other, replace)	
	
***Create Table 2. ordered logistic regression assessing associations between incentive amount and OTHERS WTA [results expressed as odds ratios]

**predictors/covariates: Language, Vaccination Status, Age, QUOTA, Gender of respondent
		
	*reorganize race var so that 'both' is the reference group
	gen race2 = 0 if QUOTA==3		//Both black and hispanic
	replace race2 = 1 if QUOTA==2	//Hispanic only
	replace race2 = 2 if QUOTA==1	//Black only

	*reorganize gender so that 'other' is the reference group 
	gen sex2 = 0 if sex ==3		//other gender
	replace sex2 = 1 if sex==1	//male
	replace sex2 = 2 if sex==2	//female


**OTHER PEOPLES' Willingness-to-Accept
*unadjusted model for OTHERS WTA
	ologit Others2 AMT, or //AIC==748.20; BIC==759.91
	margins, predict(outcome(3)) atmeans post // The probability of others being WTA = Yes given that the rest of the variables are held constant at their mean values is 47%

*adjusted models for OTHERS WTA
	*treat Age as continuous and sex as dichotomous predictor variables in the adjusted models [doing this because treating age as categorical and sex as a 3-level variable for the 2nd contigent valuation question  returned CIs that were too large to report]
	*create dichotomous sex variable
	gen female = 1 if sex ==2
	replace female = 0 if sex !=2
		
	*look at AIC values to identify model with the best fit
	ologit Others2 AMT i.Vaccine, or
	estat ic //AIC==734.42; BIC==750.04

	ologit Others2 AMT i.Vaccine Age, or
	estat ic //AIC==736.41; BIC==755.93

	ologit Others2 AMT i.Vaccine Age i.race2, or
	estat ic //AIC==739.32; BIC==766.65

	ologit Others2 AMT i.Vaccine Age i.race2 female, or
	estat ic //AIC==741.25; BIC==772.49

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or //<-- lowest AIC value, this is the final adjusted model we present in Table 2. 
	estat ic //AIC==726.86; BIC==762.01

	
*Probability that OTHERS would be WTA each amount offered in IVAC survey ($5 - $50), using adjusted model
	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=5) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=10) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=15) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=20) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=25) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=30) atmeans post

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=35) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=40) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=45) atmeans post 

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or
	margins, predict(outcome(3)) at(AMT=50) atmeans post //57.6% probability of WTA== YES

	ologit Others2 AMT, or
	margins, predict(outcome(3)) at(AMT=50) atmeans post //55.8% probability of WTA== YES
	
**What incentive amount needs to be offered per dose in order for at least 85% of OTHER PEOPLE to be WTA incentive(Others2==3) for vaccination?  
	ologit Others2 AMT, or //unadjusted model, plug and play
	margins, predict(outcome(3)) at(AMT=75 AMT=100 AMT=125 AMT=150 AMT=175 AMT=200) atmeans post //84.5% at $150 per dose and 88.7% at $175 per dose 
	marginsplot

	ologit Others2 AMT i.Vaccine Age i.race2 female language, or //adjusted model, plug and play
	margins, predict(outcome(3)) at(AMT=75 AMT=100 AMT=125 AMT=130 AMT=150 AMT=175 AMT=200) atmeans post //84.6% at $125 per dose and 85.8% at $130 per dose 
	marginsplot

	ologit Others2 AMT i.Vaccine, or //Jeff's model
	margins, predict(outcome(3)) at(AMT=75 AMT=100 AMT=125 AMT=140 AMT=150 AMT=175 AMT=200) atmeans post //85.05% of respondents WTA at $140 per dose
	marginsplot
	
label define Vaccine 0 "Unvaccinated", modify
label define Vaccine 1 "Vaccinated", modify

***Obtain 95% CIs for the amount need to get 85% of other people to be WTA incentives for vaccination 
ologit Others2 AMT i.Vaccine, or

* Save estimated coefficients in matrices
mat B = e(b)
mat W = e(V)
* Transfer to global variables b (beta), g (gamma), c (theta)
global b = B[1, 1]
global g = B[1, 3]
global c = B[1, 5]
* Get sameple mean of Vaccine (Z)
qui sum Vaccine
global Z = r(mean)
* Calculate gz (theta * Z)
global gZ = $g * $Z
* Set threshold at 85%
global p = 0.85
* Compute Y = log(p / (1-p))
global Y = log($p / (1 - $p))
* Calculate X (minimum AMT)
global X = ($Y + $c - $gZ) / $b
* Construct variance-covariance matrix V of parameters b, g, c
mat A = (1, 0, 0, 0, 0 \ 0, 0, 1, 0, 0 \ 0, 0, 0, 0, 1)
mat V = A*W*A'
* Compute first derivates dX/db, dX/dg, dX/dc
global Db = -$X / $b
global Dg = -$Z / $b
global Dc = 1 / $b
* Combine first derivates into column vector
mat D = ($Db, $Dg, $Dc)'
* Compute quadratic form variance of X
mat v = D' * V * D
* Standard error is sqrt of variance of X
global s = sqrt(v[1,1])
* Display threshold X
di $X
* Display 95% CI as +/- 1.96 * SE
di $X - 1.96 * $s, ",", $X + 1.96 * $s


**Create marginsplots for only the predictors that were significantly associated with OTHER PEOPLES' WTA: (1) vaccination status and (2) language of respondent
*Vaccination status (yes(1)/no(0))
ologit Others2 AMT i.Vaccine
margins Vaccine, predict(outcome(3)) at(AMT=(150(10)250)) vsquish //For $180 per dose, 85.4% of unvaccinated respondents and 93.4% of vaccinated respondents said others would be WTA incentives
set scheme s1mono
marginsplot, title("Predictions of Others' WTA with 95% CrIs," "by Vaccination Status of Respondent") xlabel(150(50)250) xtitle("Amount (US$)" , size(small)) ylabel(0(.5)1.0) ysc(r (0 1.0) ) ytitle("Probability Others Would Be WTA", size(small))

graph save Graph "P:\Omar Projects\Incentives to Vaccinate Against COVID-19\Final Datasets - receibved 21Sep2021\othervax.gph", replace

*Language (English(1)/Spanish(0))
ologit Others2 AMT i.language, or
margins language, predict(outcome(3)) at(AMT=(150(10)250)) vsquish //For $180 per dose, 86.1% of Spanish-speaking respondents and 92.0% of English-speaking respondents said others would be WTA incentives
set scheme s1mono
marginsplot, title("Predictions of Others' WTA with 95% CrIs," "by Primary Language of Respondent") xlabel(150(50)250) xtitle("Amount (US$)", size(small)) ylabel(0(.5)1.0) ysc(r (0 1.0) ) ytitle("Probability Others Would Be WTA", size(small))

graph save Graph "P:\Omar Projects\Incentives to Vaccinate Against COVID-19\Final Datasets - receibved 21Sep2021\otherlanguage.gph", replace


**Table 2 continued. ordered logistic regression for RESPONDENT'S Willingness-to-Accept
preserve //preserve the full dataset
	drop if You2 == . //run analysis only on the 104 unvaccinated individuals

	*unadjusted model for Respondent's WTA
	ologit You2 NEWAMT, or 
	margins, predict(outcome(3)) atmeans post // The probability of otheare held constant at their mean values is 22.7%

restore

preserve
	drop if You2 == . //run analysis only on the 104 unvaccinated individuals
	*adjusted models for Respondent's WTA
	*look at AIC values to select model with the best fit
	ologit You2 NEWAMT i.Age i.race2 i.sex2 language, or

	ologit You2 NEWAMT Age, or
	estat ic //AIC==215.5; BIC==225.09

	ologit You2 NEWAMT Age i.race2, or
	estat ic //AIC==209.12; BIC==224.99

	ologit You2 NEWAMT Age i.race2 female, or
	estat ic //AIC==209.09; BIC==227.61

	ologit You2 NEWAMT Age i.race2 female language, or //<-- final model presented in Table 2. 
	estat ic //AIC==210.71; BIC==231.87

restore

*Probability that Respondent would be WTA each amount offered in IVAC survey ($5 - $50), using adjusted models
	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =5) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =10) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =15) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =20) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =25) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =30) atmeans post

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =35) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =40) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =45) atmeans post 

	ologit You2 NEWAMT Age i.race2 female language, or
	margins, predict(outcome(3)) at(NEWAMT =50) atmeans post //25.5% probability 

	ologit You2 NEWAMT, or
	margins, predict(outcome(3)) at(NEWAMT =50) atmeans post //25.5% probability 
	
*unadjusted model
	ologit You2 NEWAMT, or
	margins, predict(outcome(3)) at(NEWAMT =50) atmeans post //29.1% probability 


**What incentive amount needs to be offered per dose in order for at least 85% of RESPONDENTS to be WTA incentive (Others2==3) for vaccination?  
	ologit You2 NEWAMT, or //unadjusted model, plug and play
	margins, predict(outcome(3)) at(NEWAMT=75 NEWAMT=100 NEWAMT=125 NEWAMT=150 NEWAMT=175 NEWAMT=200 NEWAMT=225 NEWAMT=250) atmeans post //83.0% at $200 per dose and 88.1% at $225 per dose

	ologit You2 NEWAMT Age i.race2 female language, or //adjusted model, plug and play
	margins, predict(outcome(3)) at(NEWAMT=75 NEWAMT=100 NEWAMT=125 NEWAMT=150 NEWAMT=175 NEWAMT=200 NEWAMT=225 NEWAMT=250) atmeans post //85.4% at $250 per dose 

	*calculate minimum amount needed & 95% CIs  
	ologit You2 NEWAMT 
	ereturn list
	
* Save estimated coefficients in matrices
mat B = e(b)
mat W = e(V)

* Transfer to global variables b (beta), g (gamma)
global b = B[1, 1]
global g = B[1, 3]

* Set threshold at 85%
global p = 0.85

* Compute Y = log(p / (1-p))
global Y = log($p / (1 - $p))

* Calculate X (minimum AMT)
global X = ($Y + $g) / $b

* Construct variance-covariance matrix V of parameters b, g, c
mat A = (1, 0, 0 \ 0, 0, 1 )
mat V = A*W*A'

* Compute first derivates dX/db, dX/dg, dX/dc
global Db = -$X / $b
global Dg = 1 / $b

* Combine first derivates into column vector
mat D = ($Db, $Dg)'

* Compute quadratic form variance of X
mat v = D' * V * D
* Standard error is sqrt of variance of X
global s = sqrt(v[1,1])
* Display threshold X
di $X
* Display 95% CI as +/- 1.96 * SE
di $X - 1.96 * $s, ",", $X + 1.96 * $s



***combine the two marginsplots for the stat sign predictors of OTHERS' WTA in order to create Appendix Figure 4
 gr combine "P:\Omar Projects\Incentives to Vaccinate Against COVID-19\Final Datasets - receibved 21Sep2021\othervax.gph" ///
			"P:\Omar Projects\Incentives to Vaccinate Against COVID-19\Final Datasets - receibved 21Sep2021\otherlanguage.gph" 
			**NOTE: Figure 2 was manually formatted in STATA (e.g., title, margin formatting)
			
***use Stata 16.0 to create Figure 3. Overaly two marginplots for prob WTA others and you 
	*desired model: ologit Others2 AMT
	gen pay = AMT
	gen person = Others2
	ologit person pay i.Vaccine  //Jeff's model
	margins, predict(outcome(3)) at(pay=(50(50)250)) saving(file1, replace)
	
	*desire model: ologit You2 NEWAMT
	drop pay 
	drop person
	gen pay = NEWAMT
	gen person = You2
	ologit person pay  
	margins, predict(outcome(3)) at(pay=(50(50)250)) saving(file2, replace)
	
	ssc install combomarginsplot
	set scheme s1mono
	combomarginsplot file1 file2, xtitle("Incentive Amount Per Dose ($)") label("Others" "You") noci //Titles, colors of figure manually formatted in Stata
	


			
log close
